#include <stdio.h>
#include <math.h>

#define EPSILON 1e-6  // 精度
#define MAX_ITER 1000 // 最大迭代次数

// 定义函数f(x, y)
double f(double x, double y) {
    return pow((x * x + y - 11), 2) + pow((x + y * y - 7), 2);
}

// 定义函数f(x, y)关于x和y的一阶偏导数
double df_dx(double x, double y) {
    return 4 * x * (x * x + y - 11) + 2 * (x + y * y - 7);
}

double df_dy(double x, double y) {
    return 2 * (x * x + y - 11) + 4 * y * (x + y * y - 7);
}

// 定义函数f(x, y)关于x和y的二阶偏导数
double d2f_dxx(double x, double y) {
    return 12 * x * x + 4 * y - 44;
}

double d2f_dyy(double x, double y) {
    return 4 * x + 12 * y * y - 28;
}

double d2f_dxy(double x, double y) {
    return 4 * (x + y);
}

// 牛顿法函数，用于寻找极值点
void newton_method(double *x, double *y)
{
    double hessian[2][2];
    double gradient[2];
    double hessian_inv[2][2];
    double det, dx, dy;
    int iter = 0;

    while (iter < MAX_ITER) {
        gradient[0] = df_dx(*x, *y);
        gradient[1] = df_dy(*x, *y);

        // 计算海森矩阵
        hessian[0][0] = d2f_dxx(*x, *y);
        hessian[0][1] = d2f_dxy(*x, *y);
        hessian[1][0] = d2f_dxy(*x, *y);
        hessian[1][1] = d2f_dyy(*x, *y);

        // 计算海森矩阵的行列式
        det = hessian[0][0] * hessian[1][1] - hessian[0][1] * hessian[1][0];

        // 计算海森矩阵的逆矩阵
        if (fabs(det) < EPSILON) {
            printf("行列式接近0，无法求解。\n");
            return;
        }
        hessian_inv[0][0] = hessian[1][1] / det;
        hessian_inv[0][1] = -hessian[0][1] / det;
        hessian_inv[1][0] = -hessian[1][0] / det;
        hessian_inv[1][1] = hessian[0][0] / det;

        // 更新x和y
        dx = -(hessian_inv[0][0] * gradient[0] + hessian_inv[0][1] * gradient[1]);
        dy = -(hessian_inv[1][0] * gradient[0] + hessian_inv[1][1] * gradient[1]);

        *x += dx;
        *y += dy;

        // 检查是否收敛
        if (sqrt(dx * dx + dy * dy) < EPSILON) {
            printf("收敛至极值点: x = %lf, y = %lf\n", *x, *y);
            return;
        }

        iter++;
    }

    printf("未能在最大迭代次数内收敛。\n");
}

int main()
{
    double x, y;

    // 初始化猜测值
    x = 0.0;
    y = 0.0;
    newton_method(&x, &y);

    // 可以多次尝试不同的初始值，以找到所有可能的极值点
    x = 3.0;
    y = 3.0;
    newton_method(&x, &y);

    x = -3.0;
    y = -3.0;
    newton_method(&x, &y);

    x = 3.0;
    y = -3.0;
    newton_method(&x, &y);

    x = -3.0;
    y = 3.0;
    newton_method(&x, &y);

    return 0;
}

